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Binary systems of rapidly spinning compact objects, such as black holes or neutron stars, are 
prime targets for gravitational wave astronomers. The dynamics of these systems can be very 
complicated due to spin-orbit and spin-spin couplings. Contradictory results have been presented 
as to the nature of the dynamics. Here we confirm that the dynamics - as described by the second 
post-Newtonian approximation to general relativity - is chaotic, despite claims to the contrary. 
When dissipation due to higher order radiation reaction terms are included, the chaos is dampened. 
However, the inspiral-to-plunge transition that occurs toward the end of the orbital evolution does 
retain an imprint of the chaotic behaviour. 



Gravitational wave astronomy blurs the lines between 
theory and observation by requiring accurate source mod- 
eling to facilitate detection. While it is possible to detect 
gravitational waves without precise waveform templates, 
matched filtering against a template bank is the only way 
to extract detailed information about the sources. A tem- 
plate based, matched filtering approach to gravitational 
wave data analysis is impractical if the orbital dynam- 
ics is chaotic |[ ^. Systems that exhibit sensitive 
dependence to initial conditions require template banks 
that are exponentially larger than those of non-chaotic 
systems. 

Spinning compact binaries pose a challenge to template 
based detection and parameter extraction techniques. 
The waveforms depend on a large number of parame- 
ters, including the masses of the two bodies, their spins, 
and the relative alignment of the spin and orbital angu- 
lar momentum - some 11 parameters in all. Even with a 
relatively coarse sampling of parameter space, the result- 
ing template bank can be very large. The hope is that 
hierarchical schemes can be used that start with a coarse 
sampling and proceed by successive refinement. How- 
ever, template based methods, hierarchical or otherwise, 
will not work if the underlying dynamics is chaotic. The 
sensitivity to initial conditions that characterizes chaotic 
systems ensures that waveforms that are initially nearby 
(as measured by their cross correlation over some time 
interval) will diverge exponentially with time 

A debate has arisen as to whether spinning compact 
binaries exhibit chaotic behaviour. The first indication 
of chaotic behaviour was found using a test particle ap- 
proximation [0, but chaotic orbits were only found for 
unphysically large values of the particle's spin. The prob- 
lem was also approached using the post-Newtonian ap- 
proximation to general relativity, and fractal methods 
were used to show that binaries with realistic spins ex- 
hibited chaotic behaviour 1^. Commentaries were writ- 
ten emphasizing that radiation reaction would damp the 
chaos |5[, and that the post-Newtonian approximation 
was being pushed beyond its domain of validity Q , al- 
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available. However, neither commentary disputed the 
central result of Refs. ^ - that the second post- 
Newtonian (2PN) equations of motion admit chaotic be- 
haviour. Then a paper was published "ruling out chaos 
in compact binary systems" . This study used the same 
2PN equations of motion, but a different method for es- 
tablishing chaos - Lyapunov exponents rather than frac- 
tals. The results reported in Refs.||] and sit in stark 
contrast. The trajectories that form the fractal basin 
boundaries found in belong to a set of unstable peri- 
odic orbits know as the strange repellor. These orbits 
must have positive Lyapunov exponents. Trajectories 
near the boundaries will also have positive Lyapunov ex- 
ponents, as may orbits that lie far from the boundaries. 

In what follows, we refute the claims made in Ref. 
by showing that the 2PN equations of motion do admit 
orbits with positive Lyapunov exponents. We then ex- 
plore the significance of this result by comparing three 
key timescales in the problem - the average orbital pe- 
riod To, the Lyapunov time Tx, and the decay time T^. 
If Tx is short compared to Td, the chaotic dynamics seen 
in the conservative 2PN dynamics will leave a strong im- 
print on the 2.5PN dissipative dynamics. 

The post-Newtonian equations of motion are written 
as a series expansion in v^jc?, where v is the relative 
velocity and c is the speed of light: 
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Here fj, — mim2/M is the reduced mass, M = mi -I- m2 
is the total mass, and is the relative acceleration of 
the two bodies. The product fir is given in terms of a 
series of forces, starting with the usual Newtonian force 
F^^ = 'mim2r/r^. The superscripts denote the order 



of the post-Newtonian expansion and the subscripts de- 
note the type of force. The explicit form of the higher 
order terms can be found in Refs.[^ ||, 0. Qualita- 
tively, the IPN force Fp]^ introduces perihelion preces- 
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The 2PN force Fpjy introduces isolated unsta- 
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bit (ISCO), and the possibihty of merger. The 1.5PN 
spin-orbit force F^q ' leads to precession of the orbital 
plane, as do the 2PN spin-spin F^,^ and spin- induced 

quadrupole-monopole Fgjyj forces. The spin-spin force 
is attractive for spins that are aligned and repulsive for 
spins that are anti-aligned. The 2.5PN order radiation 
reaction force, F^^\ is the first non-conservative term, 
and it causes the orbital energy E and angular momen- 
tum L to decay. Associated with the spin-orbit and spin- 
spin forces are torques that act on the spin of each body, 
causing the spins to precess - see Ref.[^ for details. While 
the expansion is known to 3PN order, we will only con- 
sider terms up to 2.5PN order to facilitate comparison 
with the results in Refs.||, |^. For the same reason, we 
also neglect the 2PN quadrupole-monopole and 2.5PN 
spin-orbit forces. In defense of these approximations, we 
point out that the terms that we keep capture the main 
qualitative features expected from full general relativity. 
If anything, the higher order non-dissipative terms are 
likely to increase the strength of the chaotic behaviour. 

The dynamics takes place in a 12-dimensional phase 
space with coordinates X — (x, p^^. Si, S2), where p^^ is 
the momentum conjugate to x, and describes the spins 
of the two bodies. In the absence of radiation reaction 
there are 6 conserved quantities, the energy E, total an- 
gular momentum J = L + Si + S2 and spin magnitudes 
Si|. Linearizing the equations of motion about a refer- 
ence trajectory X{t) gives the evolution of the difference 
SX{t) 
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The solution to this equation can be written: 
5X,{t) = U,{t)5X,{{)) . 
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The evolution matrix Lij (t) is given in terms of the linear 
stability matrix Kij by 
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with Lij{0) — Sij (repeated indices imply summation). 
The Lyapunov exponents are defined in terms of the 
eigenvalues Ai{t) of the distortion matrix A^j — LuLij: 

1 



Xi = lim ^logAi(i) . 
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The 2PN equations of motion are conservative and can be 
derived from a Hamiltonian. The expansion and vorticity 
of the flow vanishes for Hamiltonian systems (in canoni- 
cal coordinates), so that det(Ay) = 1, Aij = Aji and the 
Lyapunov exponents come in +/— pairs that measure 
the exponential shearing of the flow. The principal Lya- 
punov exponent Aj, = max(Ai) can be calculated without 
directly isolating the eigenvalues from 



= lim ^ log 
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In the limit of very long times, the principal positive 
Lyapunov exponent will dominate the trace in eqn. ^ 
By contrast, the quantity calculated in Ref.M was 



Xd = lim lim — log 

t^oodX(0)^0 t 



dX{t) 
dX{Q) 



(7) 



with dX ^ {{Xi{t) - Yi{t)){Xi{t) - Y,{t))Y/^ the Carte- 
sian distance between the 12-component vectors of a ref- 
erence trajectory X(t) and a nearby shadow trajectory 
Y{t). It must be emphasized that this is not a Lya- 
punov exponent. Eqn. ^ will automatically yield zero 
when the limit i ^ 00 is taken. However, eqn. ^ can 
represent an approximation to the Lyapunov exponent 
if an additional rescaling of the shadow trajectories is 
incorporated. The rescaling is accomplished by deter- 
mining when dX{tr) > RdX{0) for some threshold R, 
then starting a new shadow trajectory Y'(t) with initial 
conditions 



Y'{tr) = X{tr) + [Y{tr) - X {tr)) /R 



(8) 



The rescaling is repeated throughout the evolution to en- 
sure that one is accurately approximating the stability of 
the reference trajectory X{t). The problem with this 
method is that the choice of threshold can significantly 
affect the value of A^. It is possible that the apparent 
absence of rescaling in Ref. is the source of our dis- 
agreement. A far more robust method is to evolve the 
perturbation 5X(t) directly using eqn. ^j. No rescaling is 
needed as eqn. defines the dynamic stability without 
approximation. 

A second more subtle point to make regarding eqn. 
is that while dX{t) is often referred to as the "dis- 
tance between nearby trajectories in phase space", this 
statement is misleading as phase space does not admit 
a metric structure. Even if rescaled properly so that 
dX{t) w 5X{t) from eqn. the distance dX{t) only 
measures the projection of the distortion matrix onto the 
initial displacement vector: 



d^X{t) « d^{t) = dX,(0)Ay(t)dXj(0). 
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As a consequence of this additional approximation, A^ 
provides only a lower bound for Ap. 

We use three methods to estimate the principle Lya- 
punov exponent: 

Method (A) determines the full evolution matrix Lij 
from eqn. ^ and uses eqn. ^ to calculate Ap. This is 
the most numerically intensive method as it involves in- 
tegration of the 144 components of the evolution matrix 
Lij as well as the 12 components of the trajectory itself. 
The advantage of this method is that it yields an unam- 
biguous computation of the stability of an orbit with no 
approximations. 

Method (B) uses shadow trajectories and eqn. with 
a careful rescaling of the shadow orbit. This method 
involves the approximation of eqn. along with rescaling 
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FIG. 1: Fractal basin boundaries showing three possible out- 
comes for the binary system as a function of the initial spin 
alignments. 



Method (C) uses eqn. p|to evolve dX{t) along the ref- 
erence trajectory X{t) so that a total of 24 equations are 
integrated and used to calculate 
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This method combines the accuracy of integrating the 
stability equations with the approximation of projecting 
onto the distortion matrix as in eqn. ^ 

In addition to these methods we also studied the rate 
of phase decoherence in the waveforms of the reference 
and shadow trajectories. According to Ref. the phase 
difference should grow as e'^p*. In summary, all 

four methods for estimating Ap use some measure D{t), 

where D{t) is equal to ^y/{t), dX[t), d{t) or \5^{t)\, 
depending on the method. In each case, the quantity 
D(t) will have an initial power-law rise that is followed 
by exponential growth for unstable orbits. 

To illustrate the connection between the fractal struc- 
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FIG. 2: Trajectory taken from the fractal basin boundary of 
Figure 1. The axes are scaled in units of the total mass. 



ing Figure 3 of Ref. |^ in our Figure [|. The trajecto- 
ries were started in the x — y plane with initial condi- 
tions (x, x) = (S.OAf, 0, 0, 0, 0.45, 0) and spin alignments 
9i and 02 relative to the orbital angular momentum. The 
bodies have a 1 : 3 mass ratio and spins Si = 0.6mf. 
The trajectories were color coded according to their out- 
comes: Black for merger from above the x — y plane, 
dark grey for merger from below the x — y plane, white 
for more than 50 orbits, and light grey for escape be- 
yond r = lOOOAf . The lower panel in Figure |^ shows 
a detail of the fractal basin boundary, and the location 
of a long-lived orbit that lies close to the basin bound- 
ary. A portion of this trajectory is shown in Figure |^. 
The orbit has average period To = 1687Af, mean eccen- 
tricity e = 0.922 and mean semi-major axis a — 66.7Af. 
Integrating the radiation reaction force along the orbit 
gives a decay rate of < E' >= —1.26 x 10^®. In Fig- 
ure H we plot \og{D{t)/ D(G)) for this trajectory using 
methods A, B and C described above. All three methods 
yield Tx = 11500Af ~ 6.8T0 for the Lyapunov timescale, 
where we take an orbit to be a topological winding around 
the center of mass. The Lyapunov timescale is less than 
seven orbital periods, indicating that the motion is very 
chaotic. 

Similar results were found for many other orbits taken 
from Figure |l]. Most of the orbits near the boundaries 
tended to be highly eccentric (e > 0.9), by virtue of being 
on the boundary and so on the cusp between merger and 
stability. We did find some less eccentric orbits that had 
positive Lyapunov exponents. For example, the trajec- 
tory with initial conditions (x, x) = (5.5Af, 0, 0, 0, 0.4, 0), 
9i — 7r/2, 62 = 7r/6, mass ratio 1 : 3 and spins Si — ml 
is also highly chaotic. The orbit has average period 
To ~ 275A/, mean eccentricity e = 0.59 and mean semi- 
major axis a = 13.7Af. Plots of log(D(t)/i:'(0)) are 
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FIG. 3: Determining of the principle Lyapunov exponent for 
the trajectory in Figure 2. The upper Une uses method A, 
while the lower two lines (which lie over one another) use 
methods B and C. 

the phase divergence method to estimate Tx- Both meth- 
ods gave Ta = 3080M = 11.2To, which indicates that the 
orbit is highly chaotic. 

We found large numbers of orbits, with a range of mass 
ratios, spin parameters and spin alignments that had pos- 
itive Lyapunov exponents. The timescale for the chaotic 
behaviour was often a small multiple of the orbital pe- 
riod. 
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FIG. 4; Determination of the principle Lyapunov expo- 
nent for the less eccentric orbit described in the text. The 
upper line shows \og{d{t) / d{0)) while the lower line shows 
log(|5<&(t)j/|5<E.(0)|). 

To further compare with ReL §, we took up their 
case of a binary with mass ratio 1 : 1 and spins Si — mf , 
9i = 38°, 02 = 70°. We found a positive Lyapunov ex- 
ponent for (x,x) = (5.0M, 0,0, 0,0.399,0). Therefore at 
least some of the equal mass binaries demonstrate chaotic 
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jectories came with different exponents, some of which 
were zero. For example, the orbit with initial conditions 
(x,i) = (5.0M, 0,0, 0,0.428,0) gave Xp = 0. Herein lies 
an inherent weakness of the Lyapunov exponents them- 
selves. They vary from orbit to orbit. A much more 
powerful survey of the phase space scans for chaos using 
fractal basin boundaries as in Figure |^. 

We used four methods to determine Xp, along with a 
battery of numerical tests, and our results have proven ro- 
bust. We therefore confirm the chaos discovered in Refs. 

|4j , contrary to the claims of Ref . Q . It should also be 
emphasized that the fractal basin boundary method used 
m Ref. [| is an unambiguous declaration of chaos, and 
alone stands as proof of chaotic dynamics Still, 
the Lyapunov timescales can be useful for determining 
the impact of chaos on the gravitational wave detection. 

Now that we have confirmed that the 2PN dynamics 
is chaotic, we turn to the question of how significant the 
effect is. To this end we went to the next order in the 
post-Newtonian expansion and included the radiation re- 
action force. The effect of the radiation reaction force on 
the trajectory studied in Figure ^ is shown in Figure |^. 
Starting at an arbitrary point along the orbit, we see that 
the radiation reaction force drives the evolution from in- 
spiral to plunge in roughly 5 orbits. This is comparable 
to the Lyapunov timescale of roughly 11 orbits. It tells 
us is that the chaotic behaviour seen at 2PN order is 
marginal when radiation reaction is included. That is, 
chaos is damped by dissipation, but at least for some 
orbits the Lyapunov timescale is comparable to the dis- 
sipation timescale. Both the instability of the orbits and 
the degree of damping increase as merger is approached. 
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FIG. 5: A detail of trajectory (solid line) showing the effect 
of radiation reaction (dotted line) . 

There is another way to show that the chaotic be- 
haviour found in the non-dissipative 2PN dynamics does 
leave an imprint on the dissipative 2.5PN dynamics. The 
effect is illustrated in Figure |^ where trajectories with ini- 
tial conditions (x, x) — (30Af, 0, 0, 0, 0.12, 0), mass ratio 
1 : 1 and spins Si — 0.6mf were evolved for a range of 
spin-orbit alignments. The initial conditions were color 
coded using the same scheme as before. Despite the 
damping, the outcomes are intertwined in a complicated 
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system will not show fully fractal boundaries. However, 
the imprint of the underlying chaos of the conservative 
system is recorded in the amount of structure shown in 
the basin boundaries before the fractal cuts off and is 
rendered smooth. The detailed view in the lower panel 
shows that the boundaries are eventually smooth rather 
than fractal. 
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FIG. 6: Basin boundaries with radiation reaction. The upper 
panel shows a complicated intertwining of outcomes, however 
the detailed view shows that the boundaries are not fractal. 



It is worth comparing our results to the interesting 
work of Ref. Q. The authors of Ref. also found 
a positive Lyapunov exponent for spinning test parti- 
cle motion around a Schwarzschild black hole. However, 
the light companion required an unphysically large spin 
many many times maximal. Here wc find that the ad- 
ditional non-linearity from the gravitational interaction 
of the two bodies has introduced chaotic dynamics for 
physically realistic spins below maximal. We emphasize 
that the dynamics we study is only an approximation. 
We fully expect that the additional corrections at higher 
order in the PN expansion will augment the nonlinear 
behaviour and exacerbate the chaotic motion. The very 
difhculty in solving the two-body problem in general rel- 
ativity hints that the two-body problem itself, perhaps 
even without the addition of spins, is fully chaotic. 

In conclusion, there is chaos in the 2PN equations of 
motion. The chaos is damped by dissipation at 2.5PN 
order so that most orbits will only be mildly influenced 
by the complicated dynamics. What we draw from this 
is that matched filtering may survive as a viable tech- 
nique up until the innermost orbits are reached. How- 
ever, around the innermost orbits through to the plunge 
we will need to rely on other techniques (as already noted 
in Ref. ) . Importantly, the very idea of the innermost 
stable circular orbit (ISCO) must be abandoned. It is 
telling that the most unstable motion does appear to oc- 
cur in the vicinity of the homoclinic orbits. (Homoclinic 
orbits lie on the boundary between dynamical stability 
and instability The ISCO is a specific example of a 
homoclinic orbit.) The underlying chaos of the conserva- 
tive dynamics means that unstable periodic orbits crowd 
this region of phase space. The fractal basin boundaries 
are a reflection of this fractal set of unstable periodic or- 
bits. Consequently chaos can be significant for the tran- 
sition to plunge, as well as for the final orbits. 
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